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Abstract 

We study the occurrence of a Bosc-Einstcin transition in a dilute 
gas with repulsive interactions, starting from temperatures above the 
transition temperature. The formalism, based on the use of Ursell op- 
erators, allows us to evaluate the one-particle density operator with 
more flexibility than in mean-field theories, since it does not neces- 
sarily coincide with that of an ideal gas with adjustable parameters 
(chemical potential, etc.). In a first step, a simple approximation is 
used (Ursell-Dyson approximation), which allow us to recover results 
which are similar to those of the usual mean-field theories. In a second 
step, a more precise treatment of the correlations and velocity depen- 
dence of the populations in the system is elaborated. This introduces 
new physical effects, such as a marked change of the velocity profile 
just above the transition: low velocities are more populated than in an 
ideal gas. A consequence of this distortion is an increase of the crit- 
ical temperature (at constant density) of the Bose gas, in agreement 
with those of recent path integral Monte-Carlo calculations for hard 
spheres. 



1 Introduction 

The notion of Bose-Einstein condensation is not new: it was introduced 
by A. Einstein in 1925 [|J]; Bose himself played an important role in the 
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introduction of Bose-Einstein statistics, but his work was focussed on radi- 
ation (photons) - he did not generalize it to massive particles and therefore 
played no role in the discovery of the phase transition Q || . It is well known 
that, for many years, the so called "Einstein phenomenon" Q was considered 
more as a mathematical artifact of the formalism than a physical reality - 
this was for instance the case of Uhlenbeck himself |5| Q J?j] . This reaction 
is perfectly natural: indeed, the notion of accumulating particles into one 
single quantum state - with no limit on the accuracy of the definition of 
their momentum, except the size of the macroscopic container - looks rather 
paradoxical at first sight. For an ideal gas, it certainly introduces unphys- 
ical properties: density profiles which depend critically on the boundary 
conditions of the wave functions on the walls, in violation of extensivity 
||, or anomalous fluctuations of particle numbers, which become ensemble 
dependent @ @ . It is now well understood that these pathological fea- 
tures disappear as soon as some repulsive interaction between the particles 
is added; but then, by naive analogy, one could ask why the occupancy of 
a single quantum state does not disappear as well? After all, it would also 
seem perfectly natural to assume that, in an interacting system, some finite 
momentum band is highly populated^]. How can we show, from ab initio 
arguments, that the accumulation of a finite proportion of particles into a 
single quantum state is indeed a robust property against the presence of 
interactions? 

Curiously, in view of the importance of the phenomenon, the literature 
contains relatively little discussion of this question and of the possibility of 
what has sometimes be called "fractioned" or "smeared" boson condensates; 
in particular, following London's historical intuition [jl^], most textbooks 
prefer to simply assume that one single state is populated macroscopically, 
and then proceed to study the interesting consequences of this Ansatz. At 
zero temperature, a general argument was nevertheless given by Penrose and 
Onsager^ in 1956 [13]; other famous references are the work of Beliaev 14], 
of Yang [15] and of Frohlich |{l6[| ; in 1962, Girardeau 0] discussed in detail 
the possibility of what he called a "generalized condensation, where no one 
single-particle state is macroscopically occupied"; a more recent discussion 
was given by Nozieres |ll| , who concluded from a Hartree-Fock calculation at 
zero temperature that repulsive interactions tend to stabilize the single state 
occupancy. At finite temperatures, "smeared condensation" was discussed 
explicitly in an article by Luban in 1962 [16] and, more recently and in a 
mean field context, by Van den Berg et al. [2C[. There is also a large amount 

1 An elementary idea, for instance, would be to suggest that the width of the band is 
related by some inverse Fourier relation to the mean-free-path in the gas. 

2 Penrose and Onsager were the first to define the phenomenon of Bose-Einstein con- 
densation in an interacting system in terms of the eigenvalues of the one particle reduced 
density operator; in addition, from a variational argument concerning the ground state 
wave function, they give an estimation of the condensed fraction in superfluid helium four. 
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of beautiful numerical work on the subject |)21| , which leads to an impressive 
agreement with experimental data, in particular in liquid helium four; but 
one should keep in mind that numerical methods are subject to limitations 
due to finite size effects in the evaluation of narrow peaks in occupation 
numbers. For a critical discussion of experimental evidence for a condensate 
in superfluid helium four, see for instance Sokol p2| 

Two general remarks may come to mind at this stage. The first is that it 
is known, from very general considerations, that the Bose-Einstein statistics 
plays no role whatsoever at zero temperature (Boltzmann and Bose Einstein 
systems have exactly the same ground state); zero temperature arguments 
therefore do not directly address the question of how the statistics is able to 
stabilize the single state occupancy against finite temperature excitations. 
The second is that mean-field (Hartree-Fock) arguments are based on a 
approximation where the system is considered as equivalent to a gas of 
independent particles with modified energies, which automatically preserves 
the main properties of the ideal gas (for which the occurrence of Bose- 
Einstein condensation depends only on the density of states at the origin); 
they do not really address the question of the robustness of single state 
occupancy against all kinds of correlation that can be created in a system 
by interactions, but rather assume it. A theoretical study where the property 
in question would not be an ingredient, but a consequence of the results of 
the calculations, could therefore be useful. 

There are also other issues which are not completely settled, such as the 
effects of the interactions on the transition temperature itself. For instance 
one may ask if repulsive hard cores with short range will tend to increase 
or to decrease the critical temperature (at constant density]]). In 1957, the 
method of pseudopotentials [24] [25] was applied by Huang and coll. to 
the study of the properties of the phase transition in Bose hard spheres. 



3 Here, we discuss only homogeneous systems (a gas in a box), where translation in- 
variance ensures that the number density n remains constant. When this symmetry is not 
fullfilled, as is the case in atomic traps, repulsive interactions have important secondary 
effects, in particular, the spatial density of the system is changed^], with a significant 
decrease of the number density n of the atoms at the center of the trap (at constant total 
number of atoms); this reduces the degeneracy parameter nA 3 at the center of the trap and, 
obviously, the transition temperature as well. Mean-field theories can be used to account 
qualitatively for the changes of spatial distributions of atomic gases in traps and, by the 
same token, for this change of critical temperature. 

By contrast, the effects we are interested in in the present article are of a different na- 
ture; they originate from the microscopic correlations introduced between the particles by 
the iteractions, so that they remain essentally beyond the scope of mean-field approaches. 
While the latter always predict that the transition occurs when n\ 3 is equal to 2.612.. at 
the point of maximum density (center of the trap), exactly as for an ideal gas, the purpose 
of our study is precisely to study the changes of this value under the effects of interactions; 
in other words, we are interested in changes of the critical value of the degeneracy param- 
eter nA 3 , not in changes in n. Under these conditions, it is more convenient to assume 
that the gas is contained in a a box, since then n automatically remains constant, so that 
the pure correlations effects are seen with no background due to density changes. 
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The result is that the critical temperature should be slightly increased (by 
an amount proportional to the diameter a 3//2 , where a is the diameter of 
the spheres); see also ref. [26]. The physical interpretation given by the 
authors is that, by a Heisenberg type relation, "a spatial repulsion gives rise 
to a momentum space attraction", therefore facilitating the appearance of 
a condensate. But, since then, other methods of approach to the problem 
have been proposed, in particular a Hartree-Fock type approximation [pTjj. It 
turns out that this method provides the opposite result: a purely repulsive 
potential is predicted to lower the transition temperature! More precisely, 
it predicts that the change in T c depends essentially on the range of the 
potential (it vanishes for a contact potential) and that, for most of repulsive 
potentials with finite range (a "top hat potential" for instance) the exchange 
interaction tends to increase the effective mass of the particles and therefore 
to lower T c . Of course, this does not necessarily means that the two results 
are contradictory: each of them might be correct in a different domain^. 
Later, a renormalization group calculation was given by Toyoda |28[| , which 
also predicted a decrease of the critical temperature; but subsequent more 
refined calculations by Stoof and coll. led to a prediction of an increase of 
the temperature, proportional either to a |29| or to a 1 / 2 @. The most 
recent result was obtained by a numerical calculation based on the path 
integral quantum Monte-Carlo method it actually predicts the existence 
of a crossover between two regimes, at low and high densities, while at 
low densities the critical temperature is indeed increased - but even in this 
region the results do not really agree with any of the analytical calculations 
mentioned above - see also ref. [32] for another discussion of the low density 
regime. One can summarize the situation by saying that there is no present 
consensus on what is the theoretical expression of the second virial correction 
to the Bose-Einstein transition temperature in a dilute gas. 

More generally, one can be interested in a better understanding of the 
correlation that are implied by superfluidity in a gas: on a microscopic scale, 
what kind of organization in both momentum and ordinary space is responsi- 
ble for the occurrence of superfluidity? It is generally considered that super- 
fluidity is an inherently different physical phenomenon from Bose-Einstein 
condensation; but what is exactly the difference between the microscopic 
mechanisms that lead to each of these transitions, and to what extent should 
they always appear at the same time? Another related question is how the 

4 One would expect the Hartree-Fock calculation to be better for dense systems, as 
all mean-field theories; in these systems each atom interacts simultaneoulsy with many 
others, so that it can average most of its short range correlations with its neighbours, 
and experience only a mean field from them. Indeed, in a dense system such as superfluid 
helium 4, the basic prediction is the Hartree Fock calculation is borne out by experiments, 
since the critical temperature is found to be a decreasing function of pressure. 

Nevertheless the domains of validity of the two methods should overlap at least partially 
(for a potential which would be at the same time weak and with short range) so that the 
prediction of what shoud happen in this case remains unclear. 
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macroscopic wave function and its non-linear evolution |33|] |34j build up from 
individual correlated particles, in a dilute gas where atoms are free most of 
the time, collisions are binary and short, so that mean-field methods are not 
necessarily well suited? The hope is that the study of Bose-Einstein con- 
densation in dilute gases, where more precise ab initio calculations should 
be feasible than in dense systems, will allow a better understanding of the 
original phenomena observed in superfluid systems, including the stability 
against dissipation, vortices, various modes of oscillation (second sound), 
etc.; for a general review of the many recent theoretical contributions in the 
subject, see ref. [35| 

In this article, we will focus ourselves on the study of the effects of binary 
interactions on the transition temperature, in an approach of the problem 
where the populations of various individual states are kept free to vary in 
any way; they may adapt to the interactions and become more and more 
different from a Bose-Einstein distribution when the temperature of a gas is 
progressively cooled down. For this purpose, we will use the method of Ursell 
operators, already discussed in previous articles [36|, which is well adapted 
to a detailed treatment of binary short range correlations between particles 
due for instance to hard cores; for a preliminary report of the method, 
see [38]. Here, instead of calculating the partition function of the system 
as in previous work, we will find it more convenient to directly evaluate 
the one-particle density operator; it turns out that this is not only more 
direct physically, but also mathematically more convenient, mostly because 
the weights of the diagrams are much simpler than those in the partition 
function. 

It is well known [^][^][41] - and this was already emphasized in [38] in 
the context of Ursell operators - that, near the transition point, larger and 
larger exchange cycles of identical atoms become more and more important: 
infinite summations over many sizes of diagrams are therefore necessary. 
This task can be performed by using implicit integral equations; in a first 
step, we will use the simplest form for this equation, which is reminiscent of 
the Dyson equation, and this will lead us to a first form of the Ursell theory 
that is very similar to a mean- field theory. Not surprisingly then, we will 
find at this stage that the critical degeneracy parameter nA 3 is unchanged, 
keeping exactly the same value as for an ideal gas; moreover, the velocity 
distribution of atoms will also remain the same as for an ideal gas, with 
of course some change of the effective chemical potential introduced by the 
interactions. This first step is to be seen mostly as a starting base where 
known results are recovered. The second step will involve a slightly more 
elaborate integral equation for the density operator, which fortunately is not 
much more complicated than the Ursell-Dyson equation so that it can be 
solved by a similar method. The major new feature introduced is a velocity 
dependence that is no longer that of an ideal gas: it includes a distortion of 
the velocity profile, especially at low velocities, a feature which is essentially 
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beyond mean- field approximations. We will interpret this result as due to a 
spatial re-arrangement of atoms with low velocities, which allows them to 
minimize repulsion and, so to say, come closer to Bose-Einstein condensation 
than the other atoms. Consequently, the system has stronger tendency to 
populate lower states than it would have in the absence of interactions, 
which favors even more the ground state, so that it is still true that a 
single quantum state tends to be macroscopically populated; in a way, our 
reasoning can be seen as a more elaborate version of the argument concerning 
the role of interactions developed earlier by Nozieres [|l8|| , as we discuss in 
more detail below. An important consequence of this effect is that the critical 
degeneracy parameter is reduced (still assuming repulsive interactions), by 
an amount which is compatible with the numerical Monte-Carlo results of 
ref. [31]. More generally, this study provides a microscopic mechanism for 
the approach to Bose-Einstein condensation in a dilute gas. 

Section 2 below is the most technical of this article, with diagrams, 
counting, etc.; it can be skipped if the reader who is prepared to accept 
without immediate justification the integral equations proposed in section 3 
for the one-particle density operator . 



2 Ursell expansion of the one-body density oper- 
ator 

In [36] we have shown how a finite truncation of the Ursell operator series of 
the grand potential can provide virial corrections for quantum gases, even 
if they are partially degenerate; in [37] the same technique provided simple 
expressions for the one and two particle density operator. Nevertheless, in 
this article, we can not limit ourselves to these results, since they do not 
remain valid close to the Bose-Einstein condensation point. In fact any fi- 
nite truncated expression, if taken seriously, would merely exclude the phase 
transition, as discussed in § 3.3 of ref. pq ]; the transition would be replaced 
by a sharp but continuous crossover phenomenon, occurring over a finite 
range of parameters (which is independent of the size of the system but 
becomes narrower and narrower when the density of the gas decreases). But 
this conclusion arises from an incorrect simplification: it is clear that the 
validity of any truncation always ends up breaking down at some point when 
the system comes sufficiently close to Bose-Einstein condensation. Mathe- 
matically, the reason is that higher order terms in the series of perturba- 
tions, if they have smaller coefficients in a dilute system, also contain more 
denominators which diverge when the chemical potential tends to zercQ this 

5 In Ursell diagrams, every horizontal line introduces, after a summation over the size 
of the cycles, a factor (1 + /i) (where /i is the one particle distribution of the ideal gas) 
which diverges when /n 0. Therefore, the more lines the diagram contains, the more 
divergent it is. 
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implies stronger divergences near the transition point, so that higher order 
terms always become dominant at some point. Physically, the reason of 
this behavior is a divergence of the sizes of typical exchange cycles at the 
transition point [39] [40]. 

Of course, observing divergences in perturbation series at a phase transi- 
tion point is very common in physics. In the case of Bose-Einstein transition, 
nevertheless, an unusual feature is that the transition already exists in the 
ideal gas, and moreover already contains a strong singularity: the slope 
of the curve giving the density as a function of the chemical potential /i 
becomes infinite at the Bose Einstein critical point fi = (see figure 1) . 



Fig. 1: Density of an ideal g function of the chemical potential /i. If 

the effect of the interactions was simply to shift the critical value of /i by 
an amount A/i (broken lines), a first order theory would provide a 
correction which diverges when /i tends to zero. 

This makes the problem more complicated. To see why, assume for a 
moment that the effect of the interactions is just to change this critical value 
by a small amount A/i, positive if the interactions are repulsive, negative if 
they are attractive. In this case, the behavior of the density near the critical 
value of /i would be given by: 

n ~ ci - c 2 \J- (/i + A/i) (1) 

This function can be expanded in a Taylor series of powers of A/i, and the re- 
sult is a series where the term in (A/i) n is proportional to 1/ |/x| n ^ which, 
if /i tends towards zero (by negative values), diverges more and more strongly 
when its order increases; clearly the validity of the expansion breaks down 
at the transition point. This is a consequence of the non-analicity of the 
unperturbed function at the origin (an square root with an infinite deriva- 
tive) while, despite of this mathematical problem, the physical nature of the 
transition remains completely unaffected by the presence of a small correc- 
tion A/t. This example illustrates the dangers of using finite perturbation 
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series near a transition point[]; this is why we have to go beyond the expres- 
sions written in |57|] and to use a method where series are summed up to an 
infinite order. 

2.1 Operatorial derivative 

The Ursell diagrams that we will use in this article are similar to, but slightly 
different from those used in ref.[]36|]: instead of the grand potential, what will 
be expanded here is the expression of the one-particle density operator. In 
any case, the two sorts of diagrams are closely related, since reduced density 
operators can be derived from the grand canonical partition function [37] ] . 
But, instead of first introducing the diagrams and then taking derivatives, 
it turns out to be more convenient to start again the calculation from the 
beginning, mostly because this avoids introducing weights that, in a second 
step, will be cancelled in many cases. 



We use the same notation as in [36] and [37|; the P Q 's are all the N\ 



permutations of the N particles, which will be expressed as products of per- 
mutation cycles; the U^s are the Ursell operators. The first Ursell operator 
C7i is defined as function of the one particle hamiltonian Hi (kinetic energy 
for a gas in a box) as: 

E7i(l) =exp[-/3ff 1 (l)] (2) 

while the second operator U 2 is defined as a function of the two-particle 
hamiltonian H 2 (including interactions between the two particles) by: 

U 2 (l, 2) = exp [-pH 2 (l, 2)] - exp [-/3#i(l)j exp [-/?#i(2)] (3) 

(similar expressions can be written for the higher rank Ursell operators). 



Equation (9) of ref. [36] provides the following expression of the ./V particle 
partition function: 

Zn = ^J2J2 II T cluster(i,j,k...) (4) 
{Pa} {U} clusters 

where each of the ^ in equation (Q) is actually a simplified notation for two 
different summations; the sum over {U} is meant to contain a sum over all 
possible ways to write products of j7 n 's containing altogether N particles[|, 
as well as another sum over all possible non-equivalent ways to put num- 
bered particles into these U's - see equation (4) of p6|; similarly, the sum 



It is not necessarily realistic: in fact we will see in this article that the effects of 
interactions are more interesting than a simple shift of the value of the chemical potential. 
It is nevertheless interesting to note that several of the properties that we will find in our 
more elaborate study depend critically on the singularity of the curve giving the density 
of an ideal gas. 

7 This summation corresponds to the sum over the m ; 's in the notation of |36||. 
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over {P a } is meant to contain all possible ways to write products of cycles 
for N particles^ as well as all non-equivalent ways to distribute numbered 
particles in them - see §2 of [^] for more details. The T(i,j,k...) are the 
"explicit" expressions of the clusters, in fact traces involving numbered par- 
ticles (i,j,k...) which are grouped together into the same cluster by either 
Ursell operators or permutation cycles. 

To obtain the one-particle density operator, it is convenient to set (for 
n > 2): 

U n {l,2,..n) = U n (l,2,..n) x C^(l) x U x {2) x ...U^n) (5) 

where all the U n will be kept constant, while the operators U\s will be 
varied according to: 

dU 1 = dxxU 1 \tp><6\ (6) 

Here the kets | ip > and | 9 > are any kets in the one-particle state space. 
As shown in [37], the one-particle density operator p± in the canonical 



ensemble is then given by: 



(AO , Id 



< 9 I A - \«»=z- N S 2» P) 

In each term of the double sum in (Q), one has to take the derivatives of all T's 
with respect to x in succession, which amounts to taking the derivative with 
respect to the U\ operator of every numbered particle. The diagram which 
contains the particle in question then becomes an operator T c i uster (i, j, k), 
so that we obtain the expression: 

1 N 

Pi = mz ^2 Z^2^duster(hj,k)Y[T duster (m,p,q...) (8) 

N *=1 {P a } {U} rest 

where Y\ rest symbolizes the product over all remaining r c ; nster 's which have 
not be modified by the derivative; sometimes, we will call this expression 
"the triple sum". 



2.2 An example 

As an example, let us for instance take the following cluster, corresponding 
to the diagram shown in figure 2: 

r(l,2,..8) =Tri )2) ... 8 {F 2 (l,7)C/ 1 (l)C/i(2) I7i(7)I7l(8)C 6 (1 j ..6)C 2 (7, 8)} 

(9) 



This summation corresponds to the sum over the mi's in equation (7) of |36| 
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where the C"s are the permutation cycles - Cq(1, ..6) is the operator which 
creates a circular permutation of particles 1, 2, ...6 while 6*2(7, 8) is merely 
the exchange operator for particles 7 and 8. Assume for instance that we 
take the derivative with respect of particle i = 4; we then have to evaluate 
the following expression: 

rr li2j .„ 8 {F 2 (l,7)C/ 1 (l)i7 1 (2)...C/ 1 (4)C/ 1 (5)..C/ 1 (7)C/ 1 (8)C 6 (l,..6)C 2 (7,8)} 

(10) 

with the notation: 

U 1 = Ui\(p><6\ (11) 



Fig. 2: An example of an U-C term in the expansion of the partition 

function Z. 



Expression (10) is nothing but the matrix element: 

< | f duster (l,2,..8) \(p> (12) 

which can be obtained calculated by the same method of calculation as 
in p^l , with the introduction of a sufficient number of closure relations 
Y2n \ u n >< u n \ in one-particle spaces. In this way, one gets the expression: 

Em, ,n 8 < 1: %l< 2 l u n 2 I < 8 : u m | C7 2 (l, 7)C/i(l)f/i(2) x .... 

...U X {A) x C/ 1 (5)..[/ 1 (7)C/i(8) | 1 : u n2 >| 2 : u m > 
x... | 5 : u na >l 6 : u ni >| 7 : n„ 8 >| 8 : u n7 > 

(13) 

Considering first particles 1, we see that the circular permutation of indices 
in the kets at then end of the expression, together with the summation over 
indices ra 2 , 713, ...rag, introduce the product of operators: 

tr a (i,7)[z7 1 (i)] 8 tri(i)[Er 1 (i)] a (14) 

while the summation over n\ introduces a trace over particle 1; similarly, 
the summation over ras introduces the product of two additional operators 
Ui(7). Finally we get: 

Tnj {T7 2 (i, 7) [cMi)] 3 tfi(i) [^i(i)] 2 [f/i(7)] 2 } (15) 
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-which is nothing but the matrix element (the dummy index 7 is renamed 
into 2): 

<l:0|Tr 2 {[C/ 1 (l)] 2 [C/ 1 (2)] 2 I7 2 (l,2)[?7i(l)] 3 [/i(l)} | 1 : <p > (16) 

Finally we obtain: 

f ciuster (l,2,..8) =rr 2 {[C/ 1 (l)] 2 I7 2 (l,2)[[/ 1 (l)] 4 [C/ 1 (2)] 2 } 
= [C/ 1 (l)] 2 rr 2 {^ 2 (l,2)L/ 1 (2)}[[/ 1 (l)] 3 

This is, of course, an operator - no longer a number as was the initial T; but 
it can also be represented by a diagram, such as that shown in figure 3. In 
this new diagram, the lowest horizontal line no longer corresponds to a trace 
but to a product of operators ?7i(l) - as many as there are segments in this 
line (two in this particular case) - interrupted at some point by a ?7 2 (1,2), 
and then followed by another horizontal line symbolizing again the product 
of J7i's (three in this case). On the other hand, the upper horizontal line 
still corresponds to a trace over particle 2; it also contains operators U\(2), 
but here they all remain after the operator ?7 2 (1, 2). 



Fig. 3: An example of an U-C term in the expansion of the one-particle 

density-operator p 1 . 

2.3 Diagrams; weights 

More generally, all operators T in the triple sum (||) can be represented by 
diagrams where the lowest horizontal line represents a product of operators, 
beginning by either a chain of U±s, or a single U±, or actually any U n ; this 
is the main difference with diagrams contained in Z (or log Z) , where the 
initial operator was always that of highest rank. All the other horizontal 
lines correspond to traces, and begin necessarily with an operator or rank 
at least n = 2; in fact, these lines behave exactly as those in log Z. In other 
words the new diagrams are, so to say, obtained by "cutting" , or "opening" 
the old diagrams at some arbitrary point, taken in some exchange cycle[], 
which then no longer corresponds to a trace but to a product of operators - 
see the two examples shown in figure 4. 

9 This cycle may of course be of length one, corresponding to no exchange at all. 
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As in |3(J], we need a rule to decide precisely how the diagrams should be 
drawn; this is necessary in order to avoid ambiguities and double counting. 
Fortunately, the fact that the derivative has "tagged" one particle - adding 
one summation to the two which already exist in (|j) - makes the problem 
much simpler. The tagged particle determines the first operator of the lowest 
exchange cycle; it provides a well defined starting point, a "root" from which 
one can propagate horizontally along exchange cycles and vertically along 
f/2's or operators of higher rank. Actually, when only U2 operators are 
present in the diagram, no ambiguity ever occurs: the simple propagation 
from the root into the branches is sufficient to assign a well defined structure 
to each term - see for instance the first example of figure 4. Only when 
operators U n of rank n equal to 3 (or more) occur, as in the second example 
of this figure, can some ambiguity occur: which is cycle is in the middle, 
which one in the upper position? By similarity with what was done in 
[36 1, we take the following convention: beyond the first particle, which is 
determined by propagation along a previous cycle, all the other particles in 
U n 's of rank equal to 3 (or more) appear in order of increasing numbering; in 
other words, in the vertical lines, particles appear with an upward increasing 
numbering, except of course for the lowest particle which is either the tagged 
particle, or is determined from it by propagation along previous cycles and 
Un's. With this rule combined with those of |3(|, each contribution to p± 
appearing in the triple sum (||) corresponds to a perfectly well defined series 
of diagrams, starting with an operator T p , and followed by a product of 
numbers T diag , . 



Fig. 4: Two examples of diagrams occurring in the expansion of pi\ the 
diagram on the left has a weight 1, but that on the right has a weight 1/2. 

Conversely, if we did not specify any numbering in the diagrams, it is 
clear that the same series of diagrams corresponds to many different terms in 
the triple sum (pl)^; in order to get a one-to-one correspondence, we have to 
specify particle numbers at each location inside the diagrams, in other words 
to consider "numbered diagrams" where each "site" in the diagrams gets a 

10 In fact, a given T p may even originate from several different Fdiag.'s; but this is no 
longer true as soon as it contains numbered sites (the lowest of all numbers in the U„ of 
highest rank will determine the root of the T p ). 
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number. This can be done precisely by the rules mentioned in the preceding 
paragraph: starting from the tagged particle in T p , one adds horizontally the 
particles it exchanges with (in the order of the cyclic permutation) , and then 
progresses vertically along higher order rank J7 n 's: this provides a unique 
distribution of all relevant particles - as for the r^^'s, they are of course 



treated exactly as in [36]. We can then write: 



(TV) 
Pi = 



£ r P (.,.,.)II r ^.(-'-'-) ( 18 ) 



N\Z N ^ 

diag. distrib. rest 

where ^^iag stands for all the possible ways to write series of diagrams, one 
T p and one arbitrary number of r^ ag .'s, for a total number of particles equal 
to N ; in addition, ^2distrib introduces a summation over all correct ways to 
distribute numbered particles into the series. One should keep in mind that 
any random distribution of numberings in a diagram T p is not necessarily 
acceptable: in general; in fact, only a proportion g p is correct (similarly, 
only a proportion fdiag. is acceptable for diagrams T^iag., as discussed in 
||36|). We call this proportion the weight of the diagram; for instance, if a 
diagram contains one U3, the value of g p is 1/2; if it contains p operators 
t/3's its value is 1/2 P ; if it contains q operators t/4, its value is 1/6 9 , etc... 

2.4 Counting diagrams 

It is now possible to get rid of particle numbering which, of course, does not 
affect the contribution to p± of any particular term. Let us choose one given 
(non- numbered) diagram T p , call n p the number of particles it contains, and 
calculate the total coefficient that it gets from the summations of (|l8|): 



ij £ f P (.,.,.)£n r *w-(--') (19) 

distrib. diag. rest 

This can be done in two steps: first calculate the contribution of the term 
T p (i,j, k) when the particles k) are a given sub-ensemble of the ./V parti- 
cles; then sum over all the ways to select this sub-ensemble. The first step can 
be made by remarking that, for the N — n p remaining particles, the Yldistrib 
can be moved to the right of T p , which introduces the same expression as in 



J36[ ; therefore the summation of the products of T^iag. 's merely reconstructs 
the partition function Z^_ n of N — n p particles, multiplied by (N — n p )\. 
As for the n p particles, one has to take into account the number of correct 
distributions of the n p particles in the first diagram, which is (n p )\g p . For 
the second step, we first have to multiply the result by the number of ways 
to distribute n p particles among N, which is: 

N\ . , 

(20) 



(N-n p )\(n p )\ 
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Taking all these factors into account, we obtain the following result: 
1 N\ 



{n p )\ g p x(N- n p )\ Z N . np x T p (21) 



N\Z N (N - n p )\{n p )\ 
so that our final expression for the canonical ensemble is: 

p[ N) = E Zj ^9 P f P (22) 

p diagrams 

At this point, it becomes convenient to introduce the grand canonical 
ensemble and its partition function: 

Z g .c = ^2z N Z N (23) 

N 

where z = e 13 ^ is, with usual notation, the fugacity; the corresponding one- 
particle density operator p 1 (to simplify the notation, we now give up the 
index g.c): 



< | Pl | <p >= -^-^Zg.c. = ^-E^ < 9 | pf ) | <p > (24) 
We then get: 



Zg.c. dx Zg.c. N 



pi = v ^ Z Zn S 9p ^r Tp 

y N p diagrams 

= E znp 9? f P E * N - np ZN-n p 

®' c ' p diagrams N—n p 



(25) 



But the second summation reconstructs is just another expression of the 
grand canonical partition function, so that we finally obtain: 

Pi= E znp 9? f P ( 26 ) 

p diagrams 

This expression shows that the one-particle density operator is merely the 
sum of all diagrams T p , with coefficients which are the product of the weight 
g p by the fugacity raised to a power which is the number of particles con- 
tained in the diagram. 

The simplest case is, of course, the ideal gas, for which the only diagrams 
are horizontal lines; we immediately get: 

oo 

Pl = £ z n [U.T = h (27) 

71=1 
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where f\ is defined as usual by: 



Actually, even for interacting gases, this operator plays a role in the dia- 
grams, or more precisely the sum: 

oo 

J2z n [U 1 ] n = l + f 1 (29) 

?1=0 

which can be connected for instance after any operator Ui- We take the 
same convention as in |36|]: dashed horizontal lines symbolize this operator; 
the operator f\ itself is then represented by a dashed line connected to a 
one-case horizontal line (f\ = zU\ + [1 + The first terms in the Ursell 
expansion of pj are shown in figure 5; those explicitly shown in the first and 
second line correspond to the "generalized Beth-Uhlenbeck approximation 
discussed in this reference, which in terms of p 1 would lead to the expression: 

Pi= h + I* 2 [1 + /i(l)] Tr 2 {Ui(l, 2) [1 + A(2)]} [1 + h(l)] (30) 

A diagram where the same U2 re-connects to the same cycle corresponds 
to an exchange term, where U2 is replaced by its product by the exchange 
operator P ex ., as discussed in |3(|; for instance, in the second line of fig. 5, 
the second diagram is merely the exchange term of the first. 



Fig. 5: The lowest order terms in the U-C expansion of the one-particle 

density-operator p. 

3 Integral equations 

The largest eigenvalue of Ui corresponds to the ground state of the one- 
particle hamiltonian H\. When the chemical potential reaches the energy of 
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this state (zero for a gas in a box with periodic boundary conditions), the 
action of each operator zU\ leaves the ground state unchanged; then, the 
series of this operator which is contained in f\ (sum of all of its powers) 
diverges. This is not surprising since, as already mentioned, Bose-Einstein 
condensation corresponds to a situation where the size of typical exchange 
cycles is no longer limited to a few units, but increases to infinity. For an 
interacting system, the phenomenon is not limited to chains of Ux's only: 
when very long exchange cycles become important, it is expected that typical 
diagrams will also contain an arbitrary large number of operators U 2 (or 
higher rank Ursell operators) connecting these cycles. In practice, this means 
that it is not possible to limit the calculations to a few diagrams, those 
shown explicitly in figure 5 for instance; a summation over large categories 
of diagrams is indispensable. 

This is what is done in the two sections below, first within a simplified 
approach which is useful as a preliminary step (it will lead us to a mean-field 
like theory), second within a more elaborate treatment. 



3.1 A first approach: Dyson-Ursell equation 

For simplicity, we limit ourselves for the moment to terms which contain 
Ursell operators of rank 1 and 2, while a generalization to higher ranks is 
possible in a similar way. The symmetrized version of U 2 is: 

Ui = \ [1 + Pex] U 2 = \u 2 [1 + P ex ] (31) 

where P ex , is the exchange operator of the two particles contained in U 2 and 
r\ is equal to +1 for bosons (actually the only case we study here) and —1 for 
fermions; the factor 1/2 in front of this definition is introduced for the first 
part of the formula to be a projector (if necessary, one can put projectors 
on each side of XJ 2 without changing the result). 
Let us now write the following implicit equation: 

Pl (l) = /i(l) + 2z 2 [1 + Tr 2 {U 2 S (1, 2) [1 + Pl (l)] [1 + ^(2)]} (32) 

and investigate which diagrams are contained in it. The zero order term in 
U 2 is merely f±; the first order terms in U 2 are simply the second and the 



third diagram of figure 5 - in [36| we showed how terms introduced by the 
exchange operator correspond to diagrams where the same U 2 touches twice 
the same horizontal exchange cycle. More generally, it is easy to see that 

(n) 

the term p\ of order n in U 2 is obtained as a function of terms of lower 
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order by the relation^]: 

$> = 2z 2 [l + / 1 (l)]Tr 2 {^[l + / 1 (l)]^ 1) (2) 

+U!ti\l)pf- 2 \2) + + C/|pi n - 1) (l)[l + /i(2)]} 

In the right hand side, let us first consider the direct term, where is re- 
placed by Ui\ diagrammatically, in the two horizontal lines which follow C/2, 
one gets in succession all possible combinations of lower order "branches" . 
For instance, two second order terms (n = 2) are shown in the third line 
of figure 5; two others should follow with an additional U2 connected in 
two different ways to the upper cycle. As for the exchange term, it gives a 
diagram where the first U2 reconnects to the initial exchange cycle, so that 
this time it is the line between the two contact points and the line after the 
second point which play the role of the two lines after the U2; nevertheless 
the situation remains essentially the same with a different topology (taking 
again the example of second order terms, one would get 4 more diagrams 
where the first U2 is folded over the same initial cycle). We then see that, 
when the iteration is continued to infinity, it provides once, and once only, 
all possible "branched" diagrams with an arbitrary number of C/2's: indeed, 
equation ( f33|) provides a convenient way to perform a summation containing 
an infinite number of diagrams and cycles with infinite length. Nevertheless 
one should keep in mind that this is not exact either: more complicated dia- 
grams, which are not simply "branched", are still not included, for instance 
those containing loops such as those of figure 6. 



Fig. 6: Diagrams containing loops which are not included in the 
Ursell-Dyson approximation, but will be in the more elaborate integral 



equation of § 3.2; the second diagram corresponds to the exchange term of 

the first. 

Remark: equation (32) is a simple generalization of (|30[), but one could 



n The equation is written for the case where n > 1; if n = 1, the curly bracket of 
diil ) is simply equal to {U2 [1 + ,fi(l)] [1 + /i(2)]}. A more precise discussion is given in 
Appendix A 
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think of other simple iteration equations^, for instance: 

p 1 (l) = / 1 (l) + 2z 2 [l + p 1 (l)]Tr 2 {^(l,2)[l + Pl (l)][l + /0l (2)]} (34) 

(the only difference with p2] ) is that one f\ has been replaced by one pi). 
Nevertheless, one can easily convince oneself that this equation is not satis- 
factory since it just introduces redundancy, actually already at second order 
in U 2 - 



3.2 A more elaborate integral equation 

The approximation made in the preceding section can be completed by 
adding terms where U 2 operators connect twice the same pair of cycles so 
that they introduce loops. In fact, to include them, it is sufficient to add 



the following term to the right hand side of (32): 



z 



[1 + a(i)] Tr 2 { [u 2 (i, 2) [1 + Pl (i)] [1 + Pl (n] 2 [1 + Pe*: 

= 2z* [1 + Tr 2 { [U$(l, 2) [1 + Pl (l)] [1 + Pl (2)]\ 



(35) 



(the second line is easily obtained by introducing the symmetrizer [1 + P ex ] /2, 
which is a projector and can therefore be raised to any power, and using the 
commutation between this operator and U 2 ). When this is done, diagrams 
with an arbitrary number of pairs of U 2 's touching twice the same two cycles 
are included, as illustrated by the example shown in figure 7. Note that, 
since the right hand side of equation ([33]) now contains the sum of two terms, 
any combination remains possible: some cycles may be connected only once 
by a U 2 , others twice. We note in passing that, in these new diagrams, every 
U 2 does not necessarily introduce a trace over a new particle, as opposed to 
the situation in branched diagrams. 

There is no special difficulty in including now terms with an arbitrary 
number q of operators U 2 s forming a ladder between the same two cycles; 
one now has to add the following term to the right hand side of (p3|) : 



[1 + A(l)] Tr 2 { [U 2 (l, 2) [1 + Pl (l)] [1 + Pl {2)]] q [1 + P ex }} 
= 2z 2 " [1 + A(l)] Tr 2 { [Ui(l, 2) [1 + Pl (l)] [1 + Pl (2)]] q } 



(36) 



The summation of all these contributions, from q = 1 as in § |3.1| up to 
infinity, is possible and provides the following generalization of (|32|): 



p 1 (l) = / 1 (l) + 2z 2 [1 + /!(!)] Tr 2 



Uj (1,2) + [1 + Px(2)] 

-zW 2 s (1,2) [1 + Pl (l)] [1 + Pl (2)\ 

(37) 



12 Another remark is that the l's are indispensable in (|3^); otherwise diagrams where 
Us is not directly followed with U['s would be excluded. 



18 



At this level of approximation, diagrams similar to that of figure 7 with a 
completely arbitrary number of vertical lines connecting the same two cycles 
are generated. It is interesting to note that the summation of diagrams 
introduces into the trace an operator which can be seen as a two body 
generalization of the one particle distribution function of the ideal gas; one 
can easily see that more complex terms would introduce three particle terms, 
etc.. 



Fig. 7: A typical diagram containing several loops, which is included in 

equation (|37|). 

3.3 Discussion 

Clearly, the fact that almost all weights are eliminated from the diagrams 
for the density operator (if no other operator than U\ and U% is considered) 
is the source of great simplification in all calculations; similar calculations 
for the expansion of log Z would have been much more difficult. If necessary, 
one could readily include in the integral equation more complicated terms, 
such as for instance that corresponding to the diagrams of figure 8, or the 
second diagram of figure 4, chained an arbitrary number of times. We will 
see, nevertheless, that in dilute gases only a limited ensemble of diagrams 
play a role in the determination of the critical temperature. 



Fig 8: Diagrams that are not included in equation (|37|). 

4 A mean-field like approximation 

We first study equation (^2|) and calculate the value of the diagonal elements 
p k of the one-particle density operator; translational invariance in a box 
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with periodic boundary conditions ensures that the plane waves | k > are 
the eigenstates of p : . These diagonal elements are then: 



Pk =< k | Pi | k > 
while, for the ideal gas, they have the value: 

p -/3e fc 

/ k =< k I /i I k > 



1 



-Pek 



(38) 



(39) 



where f\ has been defined in (28) and is the difference between the energy 
and the chemical potential: 



2 1,2 



e* = e k - fi 



2m 



(40) 



and m the mass of the bosons. We will also replace the value of the diagonal 
matrix elements of 11% by their expression as a function of the (symmetrized) 
Ursell length afj(k), which is given in equation (17) and (4) of p3|]: 



z 2 <k,k'|[/ 2 5 |k,k' >=-^e-^ + ^' 



fe'J x afj( 



k-k 



where At is the thermal wavelength: 

At 



h 



(41) 



(42) 



(with usual notation) and V the volume of the container. But, as shown in 
this reference, the symmetric Ursell length is almost constant at low energies 
and equal to twice^] the scattering length a (if interactions occur via hard 
cores, a is merely their diameter). We therefore write: 



z 2 < k,k | U$ | k,k >- 



2^ e-^+vl x a 



(43) 



which is a very good approximation for low temperature gases (except, of 
course, if accidental collision resonances occur at very low energies); more- 
over we will assume in all the rest of this article that a is positive (repulsive 
interactions) . 



4.1 A set of non-linear coupled equations 

Projecting (|32| ) over the plane waves provides the relation: 



e -/^+vj[l + Ac ][l + /v ] (44) 



p k = /k + E[ 1 + / k ] [-^ 
k^ 

13 This factor two corresponds to the well known statistical increase of the interaction 
between two identical bosons. 
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It is convenient to introduce the variables: 



X k = e-^[1 + Pk \ (45) 
as well as the dimensionless coupling constant: 

a = 4^a (46) 

Dividing both sides of equation (^J) by the same factor [1 + f k ] then provides 
the simpler equation: 

aJ2X k 'X k -^ k X k -l = (47) 
k' 

with the notation: 

4 = 1- e^ (48) 

These equations can be obtained by stationarity conditions applied to a 
function see appendix B. 

If all variables except X k are considered as given, X k is the solution 
of a simple quadratic equation: 

a(X k ) 2 - t' k X k -1 = (49) 

with: 

C = ik~ (50) 

ancQ: 

A| fe = a £ X*> (51) 



The solution is straightforward: 



2a 



e fc + \/(4) 2 + 4a 



(52) 



but one should keep in mind that this is an implicit solution: X k depends 
on all the other X^j through A£ fe (which accounts for the mean repulsion 
exerted by all the other levels), while in turn these variables depend on X k 
through their own A£,/. The corresponding system of coupled equations 
will be solved graphically in the next section; for the moment, we assume 
for simplicity that all X^/'s are kept constant and study the variations of 



14 The bar in A£ is used to emphasize the mean-field character of this correction; in the 
next section we will include in £ fe additional terms arising from correlations. 
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Figure 9 shows this variations as a function of £ fc . We can distinguish 
between two different regimes; first, when £ fc is negative, one gets from (|52|): 



and: 



which is very reminiscent of the equations for the ideal gas; for instance, 
if (5ek <C —1, the population of the level in question is given by a Boltz- 
mann exponential. Actually, the effect of the "mean repulsion" A£ fc can be 
expressed in terms of a positive correction to the chemical potential: 

A/I fc = -/3- 1 log(l-Ae fc ) (55) 
so that, with this notation, (|54|) becomes: 



Ik = ^g.+A^ _ I = M~ &t*k) (56) 

The second regime occurs when £ fc crosses zero and becomes positive: in 
this process, the population does not diverge, in contrast to what would 
happen for the ideal gas, but reaches a new regime where its value depends 
explicitly of the coupling constant: 

t' 

X k ~ ^ (57) 
a 

Actually, since a is inversely proportional to the volume, the population 
in question becomes extensive, as expected for the ground state when it 
undergoes Bose-Einstein condensation. 



Fig. 9: variations of Xf. as a function of £j.; when £ fc is positive, reaches 
a regime where it becomes proportional to the volume of the system. 
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Equations (|5l|) and (|5^) provide the system of coupled equations to be 
solved: each depends on the value of its associated A£ fc , which in turn 
depends on all the other ^ k ''s; in general, the solution is complicated. For- 
tunately, in the thermodynamic limit, it becomes simpler since the contri- 
bution of each single X^i to A£ fc becomes negligible, so that the sum in (|5l| ) 
can be extended to include k = k and become independent of k. This is 
because, as shown by definition (|46|), the coupling constant a tends to zero 
when the volume becomes infinite, while the correction A£ fc remains finite 
(in the same limit, the discrete sum becomes an integral which contains a 
density of states which is proportional to the volume). In fact, for the con- 
tribution of a single X^ to remain significant, the level in question has to 
get a population which is proportional to the volume (extensive). 

In what follows we will therefore replace all the A£ fc 's by the same value 
A£ defined by: 

A£ = a Y, *k = |^ / rf 3 k X(k) (58) 

This is possible as long as the system is not Bose condensed; but, if one 
level ko were condensed, one would have to add its contribution aX^ to 
this integral. 



4.2 Graphical solution 

We can solve the coupled equations by a self-consistent procedure: first we 
consider A£ as given, then use its value to calculate the -Xk's, from which 
we can obtain the sum (we assume that the system is not condensed): 

^Jd^X(k)=4^J(A0 (59) 

where J is the dimensionless integral: 

Finally, we close the system by writing that the integral (59) is nothing but 
A£, which is equivalent to writing the relation: 

J(A£) = ^ A£ (61) 

and leads to a graphical solution of the self-consistent equations. 

Depending on the context, it will be convenient to express A£ either in 
terms of a change Ajl of the chemical potential: 

AJI = -/T 1 log (1 - A£) (62) 
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or of a change of the fugacity Az defined by: 



Az 



e -P*P _ i 



-zA£ (63) 



but, needless to say, these three quantities remain essentially the same pa- 
rameter. For an uncondensed level, (|53|) then provides: 

X]l ~ e f3e k + ~ I - A£ ePPk+Aft - 1 ^ 



From (|64D and ( p0|) we now get: 



J = yt^i x 9 ^ \- z ~ A ^)] (65) 



where g 3 / 2 (z) is defined as usual by: 



oo ; 

z 



33/2(z) = ^2jyt (66) 
1=1 

This function is defined only for values of z ranging from to 1. Therefore, 
when the value of \i is fixed, the integral J is defined only for values of A£ 
ranging from the minimum initial value (1 — e - ^ 4 ) up to infinity. When A£ 
takes this minimum value, the value of J goes to its maximum: 

Jm „. = £2Z2« (67) 

1 - A£ 1 - Af 

but, if one added the discrete contribution of the condensed ground state to 
the integral, the sum could take any value beyond J max . for this particular 
value of A£, the difference being proportional to the ground state population. 
Figure 10 shows plots of J as a function of A£, the chemical potential \x being 
considered as a parameter; when /i varies, according to (|67|) , the point of 
coordinates (A£, J m ax.) moves on a hyperbola (shown with a broken line). 
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Fig. 10: Plots of the function J(/x, A£) as a function of A£ for various 
values of the chemical potential \i (fi 1 < fj, 2 < fj, s < ^4); for each value of 
the chemical potential, the solution of the coupled equations is given by 
the intersection point with the straight line J = A^At/4o (point shown 
with a circle in the figure). The dotted line corresponds to the hyperbolic 
trajectory of the singular point; the critical point (// = /z 4 ) is obtained by 
intersecting this hyperbola with the same straight line (point shown with a 

square in the figure). 

For any given value of fj,, equation ( |6l| ) indicates that the value of A£ 
is obtained by intersecting the curve corresponding to J with a straight 
line going through the origin with slope Ay/4a; the intersection point is 
shown by a circle in the figure. The construction shows that, as long as is 
very negative (classical region), A£ remains only a small correction; but the 
situation changes when \x approaches zero, and more and more significant 
corrections are introduced, provided of course the dilution parameter a/ At is 
not too small. At some critical positive value JL cr , the value of A£ reaches its 
limit, and the system undergoes Bose-Einstein condensation^; this value is 
merely obtained by intersecting the hyperbola with the straight line, which 

15 The geometrical construction implies that this transition point can be reached only 
the value of the dilution parameter a/ At is sufficiently small; a simple calculation shows 
that this parameter should be smaller than [4<7 3 / 2 (l)] (which is close to 10%). If cl/\t 
was larger, the trajectory of the point with coordinates ( A£ min , J max . ) would pass above 
the straight line without ever crossing it; but a/Ar is also the parameter which determines 
the validity of the low energy approximation that we have made in the treatment of binary 
interactions, so that this situation would simply be beyond the domain of validity of the 
calculation. 



25 



provides: 

A - 4a 2.61.. 4a , , 

£cr. = T - A7 ^ 2.61..— (68) 
At 1 — A£ cr Ay 

Since at the transition point \x — A/I = 0, this corresponds to the following 
value for the chemical potential: 

10.44.. a , , 

V = tier. = A Pcr. - a T~ ( 69 ) 

p A T 



Finally, equation ( pq) shows that the number of particles as a function 
of /j, is given by: 

N = ^3 / d " k = ^3 / d " k h(k\ » ~ ( 7 °) 

At the critical value of the chemical potential, /t = A/Z, and the degeneracy 
parameter n\ T is thus predicted to have exactly the same value 2.61... as 
for the ideal gas. Therefore, what we find within this Ursell-Dyson approxi- 
mation is a relatively straightforward result: the only effect of the repulsive 
interactions is to change the effective value of the chemical potential, which 
allows positive values of /t to be reached; otherwise no change (velocity 
distribution, critical degeneracy parameter) is introduced. 



5 Beyond mean-field: velocity dependent effects 

We will now go beyond the Ursell-Dyson approximation and study which 
new effects are introduced by this more accurate treatment of the problem; 
we will see that qualitatively new effects are indeed introduced, which are 
strongly velocity-dependent. 



5.1 A new set of equations 



When (35) is added to the right hand side of (|32|), a new feature immediately 
appears: in (|35|), a momentum transfer q can now take place, introducing 
off-diagonal matrix elements^] of L^. In order to simplify the calculations, 
we will nevertheless ignore their dependence on q and merely assume that: 

, | k + q,k'-q>=-2^|e-^ 5k+ V +Sk+q+ V-J xa (71) 



z A < k,k 



V 



In appendix C we discuss these off-diagonal matrix elements and show that 
this is a reasonable approximation; but it is not essential. Equation ( |47| ) is 
then replaced by: 



a 



1 



a \ - 

2 ^ 
q 



^k+q^ k '_ f 



— £fcAk — 1 = 



(72) 



16 The matrix elements of U2 obey a momentum selection rule, but none for energy 
conservation - see for instance Appendix B. 
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while the generalization to the more complete equation 



would lead to: 



k' 



k' 



1 + 



k+q^k'-q 



-^k - ik x it -i = o 



(73) 



We note in passing that, at each "scale of the ladder" contained in the 
diagrams, momentum conservation takes place, which is why the new sum 
extends over one index q only. 

In the absence of condensation, the new term becomes an integral with 
a factor (density of states) proportional to the volume V, which cancels the 
V -1 contained in a; the result is then independent of the volume, exactly 
as in the Ursell-Dyson approximation. On the other hand, it contains an 
extra factor a/ At (dilution parameter), which ensures that the additional 
correction is smaller than that considered in the first approximation. Now, 
if the ground state is condensed, Xq becomes proportional to the volume; 
if k + k 7^ 0, the two indices k + q and k — q never vanish for the same 
value of q so that Xq now adds two constants to the integral; if k + k =0, 
the two indices vanish simultaneously when q = — k, so that one gets a 
contribution proportional to the volume. In the latter case, the sum diverges 
in the thermodynamic limit, a situation which never occurred in the Ursell- 
Dyson approximation; this is a first indication of a special coupling between 
opposite values of the momentum in the presence of condensation. 

If only one X k varies, while all the other are kept constant, the value of 
Xk is obtained as the solution of an algebraic equation; to write it explicitly, 
it is convenient to expand (|72|) term by term; when k ^ k, we obtain: 



aX k £ X k , 

kVk 



k'-q 



q^O.k^k' 



(74) 



while, if k = k, the result is changed into: 



a (X v 



1 - - (^k) 2 ~ 2Y X k+q AV 



q^O 



Adding these contributions provides the equation: 



Q 



(X k ) 4 +a(X k ) 



where £ k is now defined by: 



1 - g 12 X k+q^k- q ( X k' 
q^o k ' 



(75) 



-£ k X k -l = 
(76) 



£k = £k-aE X k'+Y£ X k' £ 4 + A- q 

k' k'^k q^0,k^k' 



(77) 
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These results generalize ( f49| ) , ( |50| ) and (^) ; they are no longer second degree 
equations, since they contain higher degree terms. Nevertheless, all these 
terms contain one extra power of a, which is inversely proportional to the 
volume, so that in the thermodynamic limit they remain negligible as long as 
Xk is not proportional (at least) to the square root of the volume - actually, 
exactly as was already the case for the quadratic term in (|49|). Under these 
conditions, the situation remains essentially similar to that discussed in the 
preceding section, and Ak is determined by the linear terms only: 

A k = ~ (78) 

The only difference is that £ k is now given by ([77]), which contains a correc- 
tion proportional to a 2 . 

If we change sums into integrals, we get: 

& = £ k -A? + <J& (79) 

where A£ is defined as above by: 

I**'**') (80) 

while the additional (second order) correction <5£ fc is given by: 

^ = 8 fe) 2 (^)V ^ x{k>) / d3 q^( k + q)^( k '-q) m 

A new feature which appears in this correction is a dependence on k, which 
will introduce velocity-dependent effects; this can be seen better by rewriting 
the second integral in the form: 

/q / ,k + k / , ,k + k / . . . 

cf 3 q X(^— +q) A(^— -q) (82) 

which appears as a convolution integral of the momentum distribution with 
itself, with a maximum when k + k = 0. As a consequence, while A k is 
coupled through 5^ k to all other A k /'s, the maximum coupling takes place 
with A_ k , which then implies that <5^ is a larger correction near the center 
of the velocity distribution (because the preferred coupling occurs with A k / 's 
which are large) than in its wings (where velocity classes are preferentially 
coupled to smaller X^j 's) . 

5.2 A shift in the transition parameter 

As in § ||, we have to solve equations consistently since the X^s depend on 
the corrections A£ and <5£j, which, in turn, depend on the A^'s, a compli- 
cated problem in general. Nevertheless, we can already get a general idea of 
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the changes introduced by the new corrections by assuming for a moment 
that they are known, or at least just the correction <5£ for the ground state. 
Figure 11 then summarizes the comparison between the Ursell-Dyson ap- 
proximation and the theory including <5£ . The full line shows the variations 
of A/Z as a function of /i, which at /x = JI cr reaches tangentially the straight 
line corresponding to the condensation condition \i + A/Z = 0; we note that, 
when fjL — > ZZ cr (by lower values), A/Z remains a function of the chemical 
potential which keeps a continuous derivative^. The broken line illustrates 
the changes to this construction introduced by <5£ , or equivalently by the 
correction 5fj, defined by analogy with ( |62[ ) (and for any value of k) as: 



The new transition point occurs when the total effective chemical potential 
of the ground state vanishes: 



which corresponds to a lower value fi cr of the chemical potential, and no 
longer to a tangential contact. An interesting feature is that the square root 
behavior of A/Z— fx near /x = ~p, cr may introduce a non-analytical behavior of 
the solution; in appendix D we write down in more detail the local expansion 
of A/Z and show, for instance, why the correction to fj, cr is second order in 
a, despite this square root behavior. Qualitatively, we can already guess 
that the critical value of the degeneracy parameter will be lowered in the 
process; this is because the excited levels get a correction <5£ fc which is smaller 
than 5^ , and therefore experience at the new transition point an effective 
chemical potential which is more negative than it was in the Ursell-Dyson 
approach (a differential effect); but we now evaluate this correction more 
qualitatively. 

17 Figure 10 shows that the slope of the curve giving the variations of J as a function of 
A£ becomes infinite when the curve connects to the vertical line - no discontinuity takes 
place in the derivative of A£ as a function of /x. 



5fi k — A/Z = /T 1 log (1 - A£ + 8£ k ) 



(83) 



li — A/x + 5/j.q = 



(84) 
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Fig. 11: Variations as a function of \x of the mean-field correction Aji to 
the chemical potential (full line); the corresponding critical value is JI cr . A 
mored detailed treatment includes the extra correction — 5fi (broken line) 
and reduces the critical value to a new value /i cr . 

A graphical solution of the consistent equations would be more compli- 
cated than in the preceding section, mostly because one single parameter 
A£ is no longer sufficient to characterize the effect of the interactions on 
all velocity classes. We can nevertheless resort to a "first <5£ iteration", 



namely calculate the <5^'s from (81) by assuming that the Xk still have the 
values that they had in the Ursell-Dyson approximation, and then use this 
correction to calculate new X^s as well as new populations. Of course, in 
theory, this should be only the initial step of an iteration procedure, which 
should be repeated until convergence is obtained, but for the moment we 
limit ourselves to this first approximation. 
Under these conditions, (|8l|) becomes: 



^ fe - 4 fe) i^) 3 ^ J ^ /i( k '^" A ^ ^(k + k';M-A7l) (85) 



where the integral J is defined by: 

3 



J(K; M -Ap) = 2(^) (^) e 2 ^ J d\ f^+^-Afl) /i(|-q;/i- 

(86) 

These two integrals can be calculated by expanding fi in series as in equation 



(27); integrating the Gaussian functions provides: 



J(K;/x-A7l) = 2 ff )e 2 ^^(/ + n)- 



3/2 x 

V ) 



l>™ n (87) 



xexp { -J^^LK 2 + p (I + n) (/i - A/Z) 
£ + n 47T 
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and: 



8£ k ^8(^-) 2 e 3 ^^[(/ + n)p + / 



nr 3/2 x 



l,n,p 

Inp A T 2 



xex M "TT^^ri - z + P(l + n + p) (ji-Aji) 
I (I + n)p + / n 47T 

The procedure is therefore the following: we start from an arbitrary 
value of fj, and use the above expressions (as well as those of § |4,2| for the 
mean-field values) to calculate the effective chemical potential ( |8"3| ) for each 
velocity class; as long as this function remains negative for all values of k, 
the population p k of each level is given by the direct generalization of fl56|): 



Pk = fifcli- A]I + 5fi k ) (89) 



Now, when \i increases, at some point condition (g4|) is met - see figure 11 - 
and condensation takes place. Figure 12 shows a plot of the variations of this 
critical value of /x as a function of the dimensionless parameter a/ Ay, with, as 
a point of comparison, the dotted straight line obtained from the mean-field 
value (|69|); we see that the mean-field theory gives a good approximation of 
the value of the chemical potential. 



Fig. 12: variations of the critical value of the chemical potential as a 
function of the dimensionless interaction parameter a/ Ay of the gas. The 

dotted upper straight line corresponds to the prediction of the 
Ursell-Dyson (mean-field) approach, equation (^); the lower line includes 
the velocity dependent effects, threated within the first iteration 
corresponding to equations (|88l) . 
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For this critical value of /i, and for three values of the parameter a/Ay, 
figure 13 then shows the variations of ~ /35[i k as a function of k, in units 
A^ 1 where the width of he thermal velocity distribution k 2 fi(k,0) is about 
4. An obvious feature in the figure is that the correction 5{A k is neither 
constant (which would lead to an unchanged critical degeneracy parameter) 
nor a simple quadratic function of k (which would lead to an effective mass 
effect). In fact it has a structure within the thermal velocity distribution, 
which resembles a Gaussian function, and which is narrower, but not by a 
large factor (2 or 3 in this case) ; this creates a intermediate physical situation 
so that quantitative predictions require a somewhat more detailed study of 
the perturbations within the velocity profile. 



Fig. 13: variations of the correction to the energy of the atoms as a 
function of their momentum, in dimensionless units kXx (the half-height 
width of the thermal distribution is then about 4); the results are obtained 
by a first iteration of the non-linear equations and the value of the chemical 
potential is adjusted to its critical value within this approximation. The 
lowest curve corresponds to cl/Xt = 0.002 and its with is 1.55; the middle 

curve to to cl/Xt = 0.003, and its width 1.75; the upper curve is for 
cl/Xt = 0.004 and its width is 1.95. This shows that the velocity-dependent 
corrections are localized at the center of the velocity profile. 

According to (|34"|), the expression of the critical density can be obtained 
by the following summation over velocities: 

n cr . = J d 3 k /i(k, n = 5fi - 6fJt k ) (90) 

This expression was used to calculate numerically the values of the critical 
degeneracy parameter nX T as a function of the dilution parameter a/Xr 
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(also equal to 1.38.. x n 1 / 3 a near the transition temperature); the results 
are shown in figure 14, expressed in terms of the difference between this 
critical value and the value 2.612.. for the ideal gas. We see that the 
integration over velocities gives rise to a correction which turns out to be, 
at least approximately, linear in a/Ay. As a point of comparison, we also 
show the results of a recent path integral quantum Monte-Carlo calculation 
[31 1; the agreement is satisfactory but not perfect; below we discuss possible 
improvements. 



Fig. 14: relative changes of the degeneracy parameter nA^ as a function of 
the density of the gas; the left vertical axis gives the opposite of this 
(negative) change, while the right vertical axis corresponds to the 
equivalent (positive) change of the critical temperature at constant density. 
The lower line shows the results of the quantum Monte Carlo calculations 
of ref. [31]; the upper curve corresponds to the results of our calculation 



within the first iteration of § 5.2; the middle curve to the more calculation 
discussed at the end of § [O], where a denominator 1 + J(k + k ) is added 
within the integral giving (but still with only one iteration of the non 

linear equations). 

At this point, it becomes possible to discuss in more detail the physical 
origin of these results. Mathematically, the key role in the mechanism is 
played by 5£ k , a correction which tends to increase effective chemical po- 
tential of a velocity class k, while A£ had just the opposite effect. In other 
words, the mean-field approximation tends to overestimate the repulsion 
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between the velocity classes, so that a correction with the opposite sign is 
introduced by (5£ k ; the correction is more important for low velocity classes 
than for rapid particles. Physically, it is indeed natural to expect that slow 
particles should be more sensitive to interactions than fast particles: they 
can more easily rearrange themselves spatially in order to minimize their 
repulsive interactions. By contrast, fast particles have too much kinetic en- 
ergy to have the same sensitivity to small perturbations; the rearrangement 
does not take place, so that they experience the full mean repulsion A£ (as 
for uncorrelated particles). The velocity dependence of the correction shown 
in figure 13; at the critical temperature, the populations of the low k levels 
are almost those of the corresponding ideal gas, with only a mass renormal- 
ization effect introduced by the curvature of the curve at the origin; but 
those of higher velocities levels get almost no correction 5£ k , so that they 
remain smaller than those of the ideal gas with chemical potential /i — Ajl. 
Integrating over all velocities provides a reduced value of the degeneracy 
parameter nA^. 

Our results go in the same direction than the views expressed by Nozieres 



[18 1 concerning the effects of repulsive interactions in a Bose gas, which 
were found to stabilize the gas against smeared condensation (fractionned 
condensate): figure 12 predicts that the populations of the low excited levels 
are indeed lower than they would be in an ideal gas, which is an indication or 
a stronger tendency to populate the ground state exclusively. The physical 
origin of the effect is nevertheless rather different: ref. (ill] is based on the 
use of a mean-field calculation, where no re-arrangement (or change of the 
velocity profile) can take place; here we find that the changes originate from 



a differential effect which is essentially beyond mean-field theory. But 18] 
treats a system which is already condensed, while we approach the transition 
from the non-condensed phase, so that the physics may well be different in 
both cases. We also note that the effective mass effect at the center of figure 
11 is of different nature than that calculated in [27]; in particular, it does 



not require the potential to have a finite range to occur (also, it has the 
opposite sign: we find a smaller effective mass than the bare mass). 

5.3 Discussion; more precise calculations 

The calculations that we have developed are not exact, for two reasons: 

(i) we have not solved the coupled equations between the <5£k's and the 
Xk, but only a "first 5t; iteration" of these equations. 

(ii) the initial equation (^|) is not exact; (|7|) is a better approximation, 
but not exact either, and can be completed with a larger class of diagrams. 



When approximation (i) is not made, equations (|85]) to (90) are replaced 



by a system containing, first, the equations providing the X^s: 

Xk " (91) 
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and, second, the equations providing the <5£ k 's (the mean repulsion A£ keeps 
the same expression as in the mean-field theory): 

tfk = 4 (^) (^) / ^ *k' >< ^( k + k ') (9 2 ) 

with: 

J(K) = 2 ^ | d 3 qI q+K/2 x X_ q+K/2 (93) 

We do not intend to solve these coupled equations in detail here, but just to 
discuss in general terms what kind of new physical effects they introduce, as 
compared to those of the preceding section. A first indication can be obtained 
by studying what would have been obtained in a "second 5^ iteration". 
Instead of the expressions ( |35"| ) and ( |36| ) containing ordinary Bose-Einstein 
distributions, we would have inside the integrals perturbed distributions, 
with already a relatively narrow enhancement near the center due to the 
first iteration. The narrowing effect of the auto-convolution integral, which 
occurred in the step, would take place again, so that we would get a still 
more localized perturbation at the center of the velocity profile. Now, when 
more and more iterations are added, narrower and narrower perturbations 
are expected to be introduced at the center, so that at the end there is no 
reason why the correction to the energy should remain quadratic at the 
origin. The situation is shown in figure 15. 



Fig. 15: when more and more iterations are performed in the calculation 
of the corrections to the energy of the particles, a sharper and sharper 
localized perturbation occurs at low velocities, so that eventually the 
spectrum no longer remains quadratic at the origin, but acquires a 

dependence in fc -3 / 2 . 
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Actually, this can also be seen from the fact that (^) gives a divergent 
integraiPl at the critical value of the chemical potential if the X^s diverge at 
the origin as k~ 2 ; such a divergence would be incompatible with the relation 
fi — Ajt + 5£ = at the transition point, as visible for instance in figure 
11. In other words, the integral giving <5£ plays the role of a restoring force 
in the problem, and the solution has to adapt to keep it finite. Assume 
then that, at the condensation point, the variations 5^ k for small values of 
k correspond to a spectrum with a power a of the momentum: 

" - *? (94) 
Then, for small velocities, the populations are proportional to: 

X k ~ k~ a (95) 
Now insert these expression into expressions (63) and (93); since 5^ k is given 



by a 6 dimension integral of three quantities which vary as k a , we get the 
relation: 

k 6~3a = k a (96) 

or: 

« = ~ (97) 

In the above reasoning, we have implicitly assumed that the behavior of 5^ k 
for small values of k depends only on the domain of integration close to the 
origin; our result is then consistent since a function in k 3 / 2 will dominate 
over a quadratic function as long as k is sufficiently small. 

As far as the second approximation is concerned, the simplest possibil- 
ity to go beyond it is to include more diagrams and to use (|73|) instead of 



(72), which amounts to adding a denominator 1 + J(k + k ) under J(k + k 



inside the integral (92). The analytical algebra then becomes even more 
complicated, but we have performed the corresponding numerical calcula- 
tions. The results are shown in the middle line in figure 14: indeed, they are 
in better agreement with the quantum Monte-Carlo results than the more 
crude approximation made before. One could therefore hope that a really 
good agreement would be attainable by including a larger class of diagrams 
in our calculations. It is also interesting to note that our conclusion concern- 
ing the exponent a = 3/2 remains valid in this case, as actually also when 
a broader class of diagrams is added in the integral equation for p 1 ; this is 
because more and more diagrams can be generated in a process where one 
U2 is added, bringing with it two X's, but also one momentum integration, 
a process in which the low k dependence of 5^ — <5£ fc is multiplied by fc 3 ~ 2a , 
a constant if a = 3/2. 

18 The first iteration (pq) of (^) already diverges, as the logarithm of \n — AjX\ |44| ; this 
divergence is taken into account in figure 10. 
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6 Conclusion 



Our study show that the effect of repulsive interactions is to shift the de- 
generacy parameter by an amount: 

x -1.3 f * -(na 3 )V3 (98) 
(n\ T )cr. At 

or, in terms of a shift of critical temperature of the gas at constant density: 

A/ ' ~ 0.7 (na 3 ) 1/3 (99) 



where a is the scattering length (assumed to be positive) of the interaction 
potential and At the thermal wavelength of the atoms. In the case of hard 
cores, these results are about 100% larger than those resulting from the nu- 



merical calculations of 1 31 ] , which provides a reasonable agreement in view 
of the approximations that we have made when solving the final non-linear 
equations; in fact, it is perfectly possible that more extensive calculations 
would lead to an even closer agreement. A striking feature of the results of 
[31 1 is the presence of two regimes, due to the increase of T c at low densities 



introduced by the repulsive interactions, followed by a decrease at higher 
densities. This is not an obvious phenomenon: usually interactions tend to 
mask quantum effects and not to enhance them, so that one could naturally 
expect a constant decrease of the transition temperature as a function of den- 
sity. Our approach explains this feature naturally, in terms of a microscopic 
re- arrangement mechanism which, in essence, automatically leads to an in- 
crease of this temperature (for repulsive interactions); the re-arrangement 
produces an enhancement of the velocity distribution at its center, which 
develops just above the transition temperature as a precursor, and allow 
the gas to reach the transition point before its density (integrated over all 
velocities) reaches the critical density of the ideal gas. In this picture, the 
mechanism for superfluidity is inherently different from Bose-Einstein con- 
densation in an ideal gas, and in a sense one could say that superfluidity 
triggers condensation at a lower density (or higher temperature) than in a 
mean field picture. At the transition point, the distortion of the velocity 
distribution bears some similarity with the results of a Bogolubov transfor- 
mation (where significant changes of the populations are also introduced for 
low energy levels), but the exponent is different, with an intermediate value 
between the high temperature and low temperature limit; in a future article, 
we intend to extend our theory to lower temperatures and condensed gases in 
order to make a more precise contact with this transformation. One can also 
remark that the singular behavior of the mean-field theory near the tran- 
sition point introduces non-analyticities in the results, mixing up together 
various orders in the natural expansion parameter a/ At; the presence of a 
square root in the calculations makes it necessary to push the calculation 
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to higher orders than one would naively expect. This suggest that theories 
using pseudopotentials may quickly reach their limit, since they are based on 
expressions of the matrix elements which are valid up to the second power 
of the potential range; from third order, the behavior of wave functions in- 
side the interaction potential (particles in the middle of a collision) becomes 
relevant, so that an Ursell approach would be safer. 

Concerning the high density regime, it remains for the moment beyond 
the domain of validity of our calculations: in denser systems, a treatment 
in terms of Ursell operators or rank exceeding 2 would become necessary. 
Doing so, one could then hope to see how a mean attractive field devel- 
ops in order to stabilize a liquid, even before Bose-Einstein condensation 
is reached - in other words, develop a more complete theory predicting the 
existence of two transitions, ordinary liquefaction and superfluid transition 
(as in helium four). But denser systems do not only differ from gases because 
of the presence of an attractive mean field; the atomic motions are also hin- 
dered by the interactions for lack of space, so that the spatial rearrangement 
which can easily take place in a gas is no longer possible. The physics of the 
effects of the interactions on exchange and on the critical temperature of 
Bose-Einstein transition is therefore significantly different. In a similar per- 
spective, we note that several of the integrals that have been introduced in 
the treatment beyond the Ursell-Dyson approximation are mathematically 
different from those of a mean-field treatment: higher dimensions integrals, 
etc.. In other words, the dimensional properties of the transition may be sig- 
nificantly altered by the interactions; this is another possibility that should 
be explored. 

Another approximation which could be released is ignoring the changes 
as a function of the relative momentum of the Ursell length (or the scattering 
length), which we have treated as exactly constant; moreover, we have not 
included the energy dependence of the off-diagonal matrix elements of XJi . A 
better treatment could for instance include a quadratic dependence on the 
relative momentum; one could then reasonably expect to recover the effective 
mass effect predicted in [p7f[ , since it is precisely due to a dependence 
of the matrix elements, and probably a better approximation to the full 
density variations predicted in [pj] . Finally, it could also be interesting 
to explore the effects of resonant changes of the cross sections at very low 
energies on the transition temperature, since they seem to be now accessible 
[45 1 [46 1|^]. One could imagine situations where sharp variations of the 



cross section could completely change the character of the Bose-Einstein 
transition, maybe even to a point where condensation would not involve a 
single level but, for instance, all levels with a given length of k; this would 
be an even more unconventional situation than the smeared condensation 
mentioned in the introduction. 
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APPENDIX A: expansion of the solution of an integral equation 



In this appendix, we study the solution of equation obtained by adding 
(36) to the right hand side of p3|): 



Pl (l) = + 2xz 2 [1 + /x(l)] Tr 2 [W(l, 2) + [W(l, 2)] 2 } (100) 

(for convenience, a variable x has been introduced in front of the trace but, 
at the end of the calculation, we can set x = 1) where: 



W(l,2) = U$(l,2)[l + Pl (l)] [1 + ^(2)] 



(101) 



We show, by a series expansion in powers of x, that this integral equation 
gives correct weights (in fact, weights equal to one) to all diagrams that 
it introduces; a similar proof could be made for the more general equation 
). If we expand the solution of ( |100D in powers of x: 



Pi 



;i) = 5>>W(i 



n=0 



(102) 



we get from fllOOD the following recurrence relation: 



pf\l) = 2z 2 [1 + h(l)] Tr 2 U 2 S (1, 2) [*,,o + #(1) 

+U 2 S "(1,2) £ 



J n—q 



+ Pi 



(2) 



i a in , 

q+q +q +q =n-l 



^,o + pS 9) (1) 



<V,o + Pi 9) ( 2 ) 



for n 7^ 0; for n = we get: 



^(1,2) 



Pi 0) (l) = /i(l) 



(103) 
(104) 
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Although equation ( |103| ) looks complicated, its content is actually simple; 
it describes how the structure of the diagrams contained in (1) is built 
from the "root" of the diagram, progressing along the lowest horizontal 
cycle/line, and adding branches made of lower orders. For example, let us 
first discuss direct terms, those for which U2 can be merely replaced by 
XJi- When the first U2 operator occurs along the lowest horizontal line, 
two cases are possible: either this operator introduces a new horizontal line 
(cycle) which is not connected anymore to the root line by any other U%, 
and this case corresponds to the first line in the right hand side of ( |103[ )); 
or it introduces a new horizontal line which is connected again by another 
U2 to the root line, corresponding to the second and third line in ( |103| ). In 
the first case, two sub-diagrams are branched after this first U2, as shown 
by the arrows in the fist diagram of fig. 16; they can occur in n possible ways 
(this corresponds to the summation over q), combining all possible lower 
orders to get a sum of orders equal to n — 1; first, the lower connection 
is replaced by a dotted line symbolizing 1 + /i(l) (the sum of ?7i(l) to all 
powers ranging from 1 to infinity) while the upper line contains the diagrams 
contained in pf" ' '; second, the lower connection gets p^ while the upper 
connection gets p™ , etc.. In the second case, the basic idea is the same, 
except that now the lower orders p± are plugged in at four different places, 
as symbolized in the second diagram of figure 16. 



Fig. 16: Diagrams summarizing the construction of p\ as a function of 

(n-l) 

Pi ■ 

As for the exchange terms, where t/.f is replaced by U^Pex., they are 
very similar, except that the first U2 reconnects to the root cycle, as shown 
in the third and fourth diagrams in figure 16. Finally, direct and exchange 
connections can be combined in any way so that, eventually, all diagrams 
corresponding to the class considered are generated. By recurrence we see 
that, if any direct Ursell diagram was contained once and only once in , 
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In) 

it will also be contained once and once only at order n in p\ . The integral 
equation therefore ascribes weights to all diagrams which are 1, as expected. 
It therefore gives an appropriate description of the approximation consid- 
ered, where two cycles are never connected more than twice. 

The generalization to any number of connections is the purpose of equa- 
tion (|37|); it can bee checked in the same way that it corresponds to appro- 
priate weights for all the diagrams. 

APPENDIX B: a potential minimization solution of the non linear 

equations 

Define the potential: 

$ = \ E a k,k'^k% - Y, (£k*k + logX k ) (105) 
k,k' k 

It is stationary under the conditions: 



which, assuming ^ 0, is equivalent to (|47D (actually to a more general 
version of these equations where all a k k / are not necessarily equal) . The 
solution of the non-linear equations of the Ursell-Dyson approximation can 
therefore be obtained by a potential minimization of the function <E>. 

Assume for instance that all Xk's except two, Xq and X\, are kept 
constant. We set: 

S = X + X 1 , , 

D = X -X 1 K ' 

and: 

f = 11 x (108) 

so that the potential function becomes, within terms which do not depend 
on Xq and Xi(from now on we assume that all a's are equal): 

$> = ^-S 2 - £S -\D- log [S 2 -D 2 ] +... (109) 

We now write: 

— — = ab — £ 



dS s S 2 -D 2 mn \ 

-A + 



dD S 2 -D 2 
The second equation provides a quadratic equation in D with the solution: 



-1± Vl + A 2 S 2 



(111) 
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but it is easy to see that only the root with the + sign in front of the radical 
is acceptable: in (f5j), the X^s are defined as positive quantities, and the 
other solution would lead to D < S. For a system contained in a box of size 
L, the energy difference between the ground state and the first excited level 
is proportional to L~ 2 , so that the same is true for A; now, if we assume that 
the sum of the populations S is macroscopic, we have S ~ L 3 and XS S> 1. 
Then: 

1 



D ~ - 



XS 1 + 



2(XS) 2 

which is, as S, proportional to L 3 . We see that: 



X 1 = (S- D)/2 
X = (S + D)/2 



L 2 
L 3 



(112) 



(113) 



which is similar to the situation for the ideal gas. Finally, S is given by the 
first equation (|110|): 



aS — £ — X 



which provides: 



S 



a 



V 



(114) 



(115) 



The advantage of the potential method is that it allows one to go beyond 
the continuous (integral) approximation of A£ and to study what happens 
to each discrete level; it would be interesting to generalize the method to 
equations (72) and ([73]) in order to check the quality of the conjectures made 
in § 5.3 concerning the linear character of the velocity profile near the center. 



APPENDIX C: matrix elements of U 2 . 

In this appendix we give a perturbative calculation of the matrix ele- 
ments of U2; for a strong potential of range 6, replacing the real potential 

the calculation remains valid up to second 



by a pseudopotential |2J] [25], 
order in b. The first order value of the part of U2 which acts in the space of 
relative motion of the two particles is: 



U r 2 eL 



13 



+ 



(116) 



where Hq is the kinetic energy hamiltonian. From this we get, still in the 
space or the relative motion: 



< x I V r 2 eL I K >= V 2 (x 



y.) — — h 



e(x') - e(x) 



(117) 
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where: 

e(x) = — (118) 

is the kinetic energy and F 2 (k) the Fourier transform of the interaction 
potential ^(r). Coming back to the full space of two interacting particles, 
we set: 

K / K / 

ki= h x ki= h x 

ko = x k 9 = x 

2 2 2 

(if x ^ x', the kinetic energy is not conserved for the relative motion, 
while it is always for the center of mass motion; in the same way, the total 
momentum is always conserved). We the have: 

flC , _ - e-/fc(*) 

< k l5 k 2 | t/ 2 (l,2) | k l5 k 2 >= e-^y 2 (q) -— + ... (120) 

e(x J - e(x) 

where is the kinetic energy of the center of mass: 

h 2 K 2 

and q is the transferred momentum: 

q= x - x ' (122) 

With this result, if we set: 

z 2 < k.k' I U„ s I k+a.k'-Q >= -4^ e ~f [*+V+*-h.H 

(123) 



z 2 < k,k I C/ 2 6 | k+q,k -q >= -4-^- e 2 L ek ^V^ k +<^V- q J x a (non-diagonal) 



we get: 



sinh /3 



e(x') — e(x) 



a (non-diagonal) = F 2 (q) , — — - — - (124) 



/3[e(x')-e(x)] 

which shows that the off diagonal matrix elements remain comparable to the 
diagonal elements within the thermal profile, with a factor of approximately 
2. 

APPENDIX D: calculating the correction to the chemical potential 
introduced by velocity-dependent effects. 

Equations (^) and ( p5| ) , together with the definition ( |62"D of Ajl, provide 
the relation: 

- 4a 1 - 

A? = ^ ^ 93/2 [z (1 - AD] (125) 
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which can also be written as: 



e -/3Ap _ e -2(3&p 



4a 
At 



53/2 



Writing the same relation for the particular value of ji: 

M t-t'cr. ^■fJ'cr. 

provides: 



e -/3A Pc _ e -2/3AjI c 



4a. 



A 



r 



where: 

Smax. = 53/2(1) = 2.612... 
Now, by difference, we obtain: 



(126) 

(127) 

(128) 
(129) 



e -/8A7» _ e -/3Ap c _ e -2/3Ap + e ~2pAjI c 



4a 
At 



53/2 
4a 



— c^P (jT) - A/J 
at 



(130) 

where the numerical coefficient c is given by (see JO] or exercise 12.3 of [24]): 

c = 3.544.. (131) 

From this we obtain: 



AtV 
4a J 



- c 



or: 



-0Aji _ e -/3Ajl c _ e -2f3Ajl + e -2/3Ajl c 



P( fl -A]l)=\^\ c- 2 [/3(A7i-A7i c )] 2 + ... 
But, to this order, we can replace A/Z by /i in the right hand side: 

/3(/i-A7J)= c- 2 tffr -A/i c )] 2 + .. • 



(132) 
(133) 

(134) 



At the critical point, fi + — A/x = so that this relation becomes: 

2 



/3o> ~ 



M \ -2 UJ .. .. \— ,]2 



4a 



or: 



4a 



/3 (Me - A*c) - T-cy/P 6 Vo 
At 



(135) 



(136) 



which is second order in o/At (since 5/U is second order in this quantity). 
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